---
title: Primary Outcome
subtitle: |
Death from any cause or requirement of new intensive respiratory
support (invasive or non-invasive ventilation) or vasopressor/inotropic
support in the 28 days after randomisation.
description: |
*This script was used to undertake the analysis of the primary outcome for ASCOT.*
author: "James Totterdell"
date: today
freeze: true
---
# Preamble
```{r}
#| label: pkgs
#| code-summary: Load packages
library (tidyverse)
library (labelled)
library (kableExtra)
library (cmdstanr)
library (posterior)
library (bayestestR)
library (bayesplot)
library (matrixStats)
library (plotly)
library (lubridate)
library (ggdist)
library (patchwork)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: analysis-sets
#| code-summary: Prepare datasets
devtools:: load_all ()
# Raw data
all_dat <- read_all_no_daily ()
all_dat_daily <- read_all_daily ()
# FAS-ITT
fas_itt_dat <- all_dat %>%
filter_fas_itt () %>%
transmute_model_cols_grp_aus_nz ()
# FAS-ITT excluding missing values
fas_itt_nona_dat <- fas_itt_dat %>%
filter (! is.na (PO))
# ACS-ITT
acs_itt_dat <- all_dat %>%
filter_acs_itt () %>%
transmute_model_cols_grp_aus_nz ()
# ACS-ITT excluding missing values
acs_itt_nona_dat <- acs_itt_dat %>%
filter (! is.na (PO))
# ACS-ITT worst case
acs_itt_wc_dat <- acs_itt_dat %>%
mutate (PO = if_else (is.na (PO) | PO == 1 , 1 , 0 ))
# ACS-ITT best case
acs_itt_bc_dat <- acs_itt_dat %>%
mutate (PO = if_else (is.na (PO) | PO == 0 , 0 , 1 ))
# ACS-PP - exclude participants not per-protocol
pp <- read_per_protocol_list ()
all_pp <- all_dat %>% left_join (pp, by = "StudyPatientID" )
acs_pp_dat <- all_pp %>%
filter_acs_itt () %>%
filter (is.na (Reason)) %>%
transmute_model_cols_grp_aus_nz ()
acs_pp_nona_dat <- acs_pp_dat %>%
filter (! is.na (PO))
# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
filter_concurrent_intermediate ()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>%
filter (! is.na (PO))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
filter_concurrent_std_aspirin ()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>%
filter (! is.na (PO))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
filter_concurrent_therapeutic ()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>%
filter (! is.na (PO))
```
```{r}
#| label: stan-models
#| code-summary: Load models
logistic_mod <- cmdstan_model ("stan/binary/logistic.stan" )
logistic_site_mod <- cmdstan_model ("stan/binary/logistic_site.stan" )
logistic_epoch_mod <- cmdstan_model ("stan/binary/logistic_epoch.stan" )
logistic_site_epoch_llik <- cmdstan_model ("stan/binary/logistic_site_epoch_loglik.stan" )
logistic_site_epoch_llik_int <- cmdstan_model ("stan/subgroup/logistic_site_epoch_loglik_shrinkage.stan" )
primary_mod <- cmdstan_model ("stan/binary/logistic_site_epoch.stan" )
```
```{r}
#| label: helper-functions
#| code-summary: Helper functions
make_stan_data <- function (
dat,
vars,
outcome,
beta_sd = c (2.5 , rep (1 , nXassign), 2.5 , 10 ),
ctr = contr.orthonorm)
{
X <- model.matrix (
as.formula (paste ("~" , paste (vars, collapse = " + " ))),
data = dat,
contrast = list (
AAssignment = ctr,
CAssignment = ctr)
)
nXassign <- sum (grepl ("Assignment" , colnames (X)))
y <- dat[[outcome]]
N <- dim (X)[1 ]
K <- dim (X)[2 ]
epoch <- dat$ epoch
M_epoch <- max (dat$ epoch)
region <- dat[["ctry_num" ]]
M_region <- max (region)
site <- dat[["site_num" ]]
M_site <- max (site)
region_by_site <- region_by_site <- dat %>%
dplyr:: count (ctry_num, site_num) %>%
pull (ctry_num)
list (N = N, K = K, X = X, y = y,
M_region = M_region, region = region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
beta_sd = beta_sd)
}
# Always includes the assigned interventions
make_primary_model_data <- function (
dat,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
ctr = contr.orthonorm) {
# Domain specific intercept as in interims
XA <- make_domA_design (dat, ctr)
# Domain specific intercept as in interims
XC <- make_domC_design (dat, ctr)
if (any (XC[, 1 ] != 1 )) {
add_intercept = TRUE
} else {
add_intercept = FALSE
}
if (! is.null (vars)) {
Xother <- model.matrix (
as.formula (paste (" ~ " , paste (vars, collapse = " + " ))),
data = dat
) [, - 1 ]
} else {
Xother <- NULL
}
X <- cbind (XC, XA, Xother)
if (add_intercept) {
X <- cbind (randPlatform = 1 , X)
}
nXtrt <- sum (grepl ("rand" , colnames (X))) - 1
epoch <- dat[["epoch" ]]
M_epoch <- max (epoch)
region <- dat[["ctry_num" ]]
M_region <- max (region)
site <- dat[["site_num" ]]
M_site <- max (site)
region_by_site <- region_by_site <- dat %>%
dplyr:: count (ctry_num, site_num) %>%
pull (ctry_num)
y <- dat[["PO" ]]
N <- dim (X)[1 ]
K <- dim (X)[2 ]
beta_sd <- c (2.5 , rep (1 , nXtrt), beta_sd_var)
out <- list (
N = N, K = K, X = X, y = y,
M_region = M_region, region = region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
beta_sd = beta_sd,
ctrC = attr (XC, "contrasts" )[[1 ]],
ctrA = attr (XA, "contrasts" )[[1 ]]
)
return (out)
}
fit_primary_model <- function (
dat,
model = primary_mod,
vars = NULL ,
beta_sd_var = NULL ,
ctr = contr.orthonorm,
seed = 32915 ,
...
) {
mdat <- make_primary_model_data (
dat,
vars,
beta_sd_var,
ctr
)
snk <- capture.output (
mfit <- model[["sample" ]](
data = mdat,
seed = seed,
refresh = 0 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
chains = 8 ,
parallel_chains = min (8 , parallel:: detectCores ()),
...
))
mpars <- mfit$ metadata ()$ model_params
keep <- mpars[! grepl ("(_raw|epsilon_)" , mpars)]
mdrws <- as_draws_rvars (mfit$ draws (keep))
names (mdrws$ beta) <- colnames (mdat$ X)
if (any (grepl ("site" , keep))) {
site_map <- dat %>% dplyr:: count (site_num, site)
names (mdrws$ gamma_site) <- site_map$ site
}
if (any (grepl ("epoch" , keep))) {
epoch_map <- dat %>% dplyr:: count (epoch, epoch_lab)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
}
# Transformed samples
Ca <- mdat$ ctrA
Cc <- mdat$ ctrC
# Get constrained parameters from contrast transformation
mdrws$ Acon <- Ca %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- as.vector (Cc %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))])
names (mdrws$ Ccon) <- rownames (Cc)
# Get treatment effect in case of non treatment-coding
mdrws$ Atrt <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
# Get odds ratios
mdrws$ OR <- exp (mdrws$ Ctrt)
mdrws$ OR <- c (mdrws$ OR, exp (mdrws$ beta[! grepl ("(Intercept|rand)" , names (mdrws$ beta))]))
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
make_po_table <- function (dat, format = "html" ) {
tab <- dat %>%
mutate (CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C3" , "C4" ),
labels = intervention_labels ()$ CAssignment[- 1 ])
) %>%
group_by (CAssignment) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (PO)), 100 * mean (is.na (PO))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (PO)), 100 * mean (! is.na (PO))),
` Met primary outcome ` = sprintf ("%i (%.1f)" , sum (PO, na.rm = TRUE ), 100 * mean (PO, na.rm = TRUE )),
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Randomised = n (),
` Outcome missing ` = sprintf ("%i (%.1f)" , sum (is.na (PO)), 100 * mean (is.na (PO))),
` Outcome observed ` = sprintf ("%i (%.1f)" , sum (! is.na (PO)), 100 * mean (! is.na (PO))),
` Met primary outcome ` = sprintf ("%i (%.1f)" , sum (PO, na.rm = TRUE ), 100 * mean (PO, na.rm = TRUE )),
)
) %>%
mutate (CAssignment = fct_inorder (CAssignment)) %>%
gather (key, value, - CAssignment, factor_key = T) %>%
spread (CAssignment, value)
colnames (tab)[1 ] <- "n ( \\ %)"
if (format == "latex" ) {
colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )
}
kable (tab,
format = format,
align = "lrrrrr" ,
escape = FALSE ,
booktabs = TRUE
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position" )
}
summarise_posterior <- function (dat, futval = 1.1 ) {
dat %>%
mutate (
Median =
sprintf ("%.2f" , median (Posterior)),
` 95 \\ % CrI ` =
sprintf ("(%.2f, %.2f)" , quantile (Posterior, 0.025 ), quantile (Posterior, 0.975 )),
` Mean (SD) ` =
sprintf ("%.2f (%.2f)" , mean (Posterior), sd (Posterior)),
` Pr(OR < 1) ` =
sprintf ("%.2f" , Pr (Posterior < 1 )),
` Pr(OR > 1/1.1) ` =
sprintf ("%.2f" , Pr (Posterior > 1 / futval))
)
}
make_subgroup_summary_table <- function (dat, format = "html" ) {
L <- nrow (dat)
row_st <- seq (1 , L, by = L / 3 )
row_ed <- seq (L/ 3 , L, by = L / 3 )
kable (
dat[, - 1 ] %>% select (- Posterior),
format = format,
booktabs = TRUE ,
align = "lrrrrr" ,
escape = FALSE
) %>%
kable_styling (
latex_options = "HOLD_position" ,
font_size = 9
) %>%
group_rows ("Intermediate" , row_st[1 ], row_ed[1 ]) %>%
pack_rows ("Low with aspirin" , row_st[2 ], row_ed[2 ]) %>%
pack_rows ("Therapeutic" , row_st[3 ], row_ed[3 ])
}
plot_or_densities <- function (rvs) {
tibble (Contrast = fct_inorder (names (rvs)), RV = rvs) %>%
ggplot (., aes (y = Contrast, xdist = RV)) +
stat_halfeye (
aes (fill =
stat (cut_cdf_qi (
cdf,
.width = c (.5 , .8 , .95 , 0.99 ),
labels = scales:: percent_format ()))),
adjust = 1 , n = 1001 , .width = c (0.5 , 0.8 , 0.95 )
) +
scale_fill_brewer (
palette = "Reds" ,
direction = - 1 ,
na.translate = FALSE ) +
labs (
x = "Odds ratio contrast" ,
fill = "Interval"
) +
scale_x_log10 (breaks = c (0.1 , 0.25 , 0.5 , 1 , 2 , 4 , 10 )) +
geom_vline (xintercept = 1 )
}
odds_ratio_summary_table <- function (OR, format = "html" , fn = NULL ) {
out <- tibble (
Parameter = names (OR),
Median = median (OR),
` 95% CrI ` = apply (
quantile (OR, c (0.025 , 0.975 )), 2 ,
function (x) sprintf ("(%.2f, %.2f)" , x[1 ], x[2 ])),
` Mean (SD) ` = sprintf ("%.2f (%.2f)" , E (OR), sd (OR)),
` Pr(OR < 1) ` = Pr (OR < 1 ),
) %>%
kable (
format = format,
digits = 2 ,
align = "lrrrr" ,
linesep = "" ,
booktabs = TRUE ) %>%
kable_styling (
font_size = 9 ,
bootstrap_options = "striped" ,
latex_options = "HOLD_position" )
if (! is.null (fn) & format == "latex" ) {
save_tex_table (out, fn)
} else {
return (out)
}
}
decision_quantity_summary_table <- function (OR, format = "html" , fut_val = 1.1 ) {
tdat <- tibble (
Intervention = names (OR),
posterior = OR,
Posterior = sprintf ("%.2f (%.2f, %.2f)" ,
median (posterior),
quantile (posterior, 0.025 ),
quantile (posterior, 0.975 )),
` Superior \n Pr(OR = min(OR)) ` = sprintf ("%.2f" , E (posterior == rvar_min (posterior))),
` Effective \n Pr(OR < 1) ` = sprintf ("%.2f" , Pr (posterior < 1 )),
` Futile \n Pr(OR > 1/1.1) ` = sprintf ("%.2f" , Pr (posterior > 1 / fut_val)),
` Equivalent \n Pr(1/1.1 < OR < 1.1) ` = sprintf ("%.2f" , Pr (posterior < fut_val & posterior > 1 / fut_val))
) %>%
select (- posterior)
tdat[1 , 2 ] <- "1.00"
tdat[1 , - (1 : 3 )] <- "-"
if (format == "latex" ) {
colnames (tdat) <- linebreak (colnames (tdat), align = "c" , linebreaker = " \n " )
}
out <- kable (
tdat,
format = format,
digits = 2 ,
align = "lrrrrr" ,
linesep = "" ,
escape = F,
booktabs = TRUE ) %>%
kable_styling (
font_size = 9 ,
bootstrap_options = "striped" ,
latex_options = "HOLD_position"
)
return (out)
}
```
# Outline
This document presents the results of analyses of the primary outcome. The first section introduces the primary outcome definition and it's derivation from the database fields. Descriptive analyses are then presented followed by the model-based analyses.
For the model-based analyses, results are presented for each population and sub-model.
# Primary Outcome Definition
The primary outcome is a composite of death, need for new respiratory support, or vasopressor/inotropic support. From the protocol:
> Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation. This includes any participant who receives non-invasive mechanical ventilation (either CPAP or BIPAP, apart from the below considerations) any time after enrolment even if not transferred to ICU. It does NOT include the use of humidified high-flow nasal prong oxygen.
>
> Participants on pre-existing home BiPAP or CPAP will not be considered to have met the primary outcome unless they have either:
>
> - required invasive mechanical ventilation (i.e. intubation), or
> - graduated from CPAP only whilst asleep to BiPAP at any time, or
> - graduated from BiPAP only whilst asleep to BiPAP for \> 12 hours/day, or
> - died by day 28
>
> There may be cases where a patient has been assessed as requiring intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support, but the patient or family declined treatment and the patient was discharged home. If attempts to obtain 28-day data are unsuccessful or not possible, and the investigator had deemed at the time of discharge that the patient would be highly likely to die within 28 days from randomisation, these participants will be deemed to have met the primary outcome.
# Derivation of the Outcome
Derivation of the outcome requires checking of the daily, discharge, and day 28 data extracts. On the daily data, there is a variable `DD_PrimaryEndpointReachedToday` however this was coded incorrectly in the original database and therefore fails to capture some participants. Additionally, given the composite nature of the outcome it is useful to check all components individually as well as the composite outcome. Therefore, this variable is not used in the derivation, but is included as a cross-check.
Below, each component is summarised in aggregate.
## Day 28 mortality
```{r}
#| label: po-component-death
#| code-summary: Summarise the death component
mort_table <- all_dat %>%
filter_acs_itt () %>%
group_by (CAssignment) %>%
summarise (
` Mortality_Alive at day 28_ ` = sum (1 - D28_death, na.rm = TRUE ),
` Mortality_Death within 28 days_Total ` = sum (D28_death, na.rm = TRUE ),
` Mortality_Death within 28 days_Prior to discharge ` = sum (DIS_death & DIS_day <= 28 , na.rm = TRUE ),
` Mortality_Death within 28 days_Post-discharge ` = sum (D28_death == 1 & DIS_death == 0 , na.rm = TRUE ),
` Mortality_Unknown_ ` = sum (is.na (D28_death))
) %>%
gather (key, Frequency, - CAssignment, factor_key = T) %>%
spread (CAssignment, Frequency) %>%
separate (key, into = c ("Measure" , "Outcome" , "Breakdown" ), sep = "_" ) %>%
mutate (
Overall = C1 + C2 + C3 + C4,
across (C1: Overall, ~ sprintf ("%i (%.0f)" , .x, 100 * .x / sum (.x[c (1 ,2 ,5 )])))
)
kable (
mort_table %>% select (- 1 ),
caption = "Composite primary endpoint breakdown - mortality." ,
align = "llrrrrr" ) %>%
kable_styling (font_size = 11 ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'custom' , custom_latex_hline = 1 ) %>%
group_rows ("Mortality" , 1 , nrow (mort_table))
```
## Vasopressor/inotropic requirement
```{r}
#| label: po-component-vasopressor
#| code-summary: Summarise the vasopressor/inotropes component
vaso_table <- all_dat %>%
filter_acs_itt () %>%
group_by (CAssignment) %>%
summarise (
` Vasopressor/inotropes_Not required_ ` = sum (1 - ANY_vasop, na.rm = TRUE ),
` Vasopressor/inotropes_Use within 28 days_Total ` = sum (ANY_vasop, na.rm = TRUE ),
` Vasopressor/inotropes_Use within 28 days_Prior to discharge ` = sum (DD_vasop, na.rm = TRUE ),
` Vasopressor/inotropes_Use within 28 days_Post-discharge ` = sum (D28_vasop, na.rm = TRUE ),
` Vasopressor/inotropes_Unknown_ ` = sum (is.na (ANY_vasop))
) %>%
gather (key, Frequency, - CAssignment, factor_key = T) %>%
spread (CAssignment, Frequency) %>%
separate (key, into = c ("Measure" , "Outcome" , "Breakdown" ), sep = "_" ) %>%
mutate (
Overall = C1 + C2 + C3 + C4,
across (C1: Overall, ~ sprintf ("%i (%.0f)" , .x, 100 * .x / sum (.x[c (1 ,2 ,5 )])))
)
kable (
vaso_table %>% select (- 1 ),
caption = "Composite primary endpoint breakdown - vasopressor/inotropes." ,
align = "llrrrr" ) %>%
kable_styling (font_size = 11 ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'custom' , custom_latex_hline = 1 ) %>%
group_rows ("Vasopressor/inotropes" , 1 , nrow (vaso_table))
```
## New intensive respiratory support
```{r}
#| label: po-component-respiratory
#| code-summary: Summarise the ventilation component
vent_table <- all_dat %>%
filter_acs_itt () %>%
group_by (CAssignment) %>%
summarise (
` Ventilation_Not required_ ` = sum (1 - ANY_vent, na.rm = TRUE ),
` Ventilation_Use within 28 days_Total ` = sum (ANY_vent, na.rm = TRUE ),
` Ventilation_Use within 28 days_Prior to discharge ` = sum (ANY_DD_vent, na.rm = TRUE ),
` Ventilation_Use within 28 days_Post-discharge ` = sum (ANY_vent & ANY_DD_vent == 0 , na.rm = TRUE ),
` Ventilation_Unknown_ ` = sum (is.na (ANY_vent))
) %>%
gather (key, Frequency, - CAssignment, factor_key = T) %>%
spread (CAssignment, Frequency) %>%
separate (key, into = c ("Measure" , "Outcome" , "Breakdown" ), sep = "_" ) %>%
mutate (
Overall = C1 + C2 + C3 + C4,
across (C1: Overall, ~ sprintf ("%i (%.0f)" , .x, 100 * .x / sum (.x[c (1 , 2 , 5 )])))
)
kable (
vent_table %>% select (- 1 ),
caption = "Composite primary endpoint breakdown - respiratory support." ,
align = "llrr" ) %>%
kable_styling (font_size = 11 ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'custom' , custom_latex_hline = 1 ) %>%
group_rows ("Ventilation" , 1 , nrow (vent_table))
```
## DAMA and Likely to Die
There were 3 participants who were DAMA and identified as likely to die. However, the day 28 status was observed for all 3 participants, therefore, no-one met the primary endpoint due to DAMA and unknown day 28 status.
```{r}
#| label: po-components-dama
#| code-summary: Summarise the DAMA component
all_dat %>%
filter_acs_itt () %>%
select (DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
dplyr:: count (DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
spread (D28_death, n) %>%
kable (
col.names = c ("DAMA" , "Likely to die by day 28" , "No" , "Yes" , "(Missing)" )
) %>%
kable_styling () %>%
add_header_above (c (" " = 2 , "Death by day 28" = 3 ))
```
## Composite
```{r}
#| label: tbl-po-composite
#| code-summary: Summarise the composite
#| tbl-cap: |
#| Composite primary endpoint breakdown - overall.
comp_table <- all_dat %>%
filter_acs_itt () %>%
group_by (CAssignment) %>%
summarise (
` Primary Outcome_No_ ` = sum (1 - PO, na.rm = TRUE ),
` Primary Outcome_Yes_ ` = sum (PO, na.rm = TRUE ),
` Primary Outcome_Unknown_Total ` = sum (is.na (PO)),
` Primary Outcome_Unknown_Day 28 status ` = sum (is.na (PO) & is.na (D28_death)),
` Primary Outcome_Unknown_Vasopressor/inotropes ` =
sum (is.na (PO) & ! is.na (D28_death) & is.na (ANY_vasop))
) %>%
gather (key, Frequency, - CAssignment, factor_key = T) %>%
spread (CAssignment, Frequency) %>%
separate (key, into = c ("Measure" , "Outcome" , "Breakdown" ), sep = "_" ) %>%
mutate (
Overall = C1 + C2 + C3 + C4,
across (C1: Overall, ~ sprintf ("%i (%.0f)" , .x, 100 * .x / sum (.x[1 : 3 ])))
)
kable (
comp_table %>% select (- 1 ),
align = "llrr" ) %>%
kable_styling (font_size = 11 ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'custom' , custom_latex_hline = 1 ) %>%
group_rows ("Primary outcome" , 1 , nrow (comp_table))
```
```{r}
#| label: tbl-po-composite-all
#| code-summary: Summarise the all components
#| tbl-cap: |
#| Composite primary endpoint breakdown.
tab <- bind_rows (comp_table, mort_table, vaso_table, vent_table)
tab_head <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (CAssignment = factor (CAssignment, labels = intervention_labels ()$ CAssignment[- 1 ])) %>%
mutate (label = paste0 (CAssignment, "<br>(n = " , n, ")" )) %>%
pull (label) %>%
c (paste0 ("Overall<br><br>(n = " , nrow (all_dat %>% filter_acs_itt ()), ")" ))
colnames (tab)[4 : 8 ] <- tab_head
kable (
tab %>% select (- 1 ),
align = "llrrrrr" ,
escape = F) %>%
kable_styling (font_size = 11 ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'custom' , custom_latex_hline = 1 ) %>%
group_rows ("Primary outcome" , 1 , nrow (comp_table)) %>%
group_rows ("Mortality" , nrow (comp_table) + 1 , nrow (comp_table) + nrow (mort_table)) %>%
group_rows ("Vasopressor/inotropes" , 11 , 15 ) %>%
group_rows ("Ventilation" , 16 , 20 )
```
```{r}
#| label: save-composite-primary-table
#| code-summary: Save table to outputs
colnames (tab)[4 : 8 ] <- linebreak (tab_head, align = "c" , linebreaker = "<br>" )
out <- kable (
tab %>% select (- 1 ),
format = "latex" ,
align = "llrrrrr" ,
booktabs = TRUE ,
escape = F) %>%
kable_styling (font_size = 9 , latex_options = "HOLD_position" ) %>%
collapse_rows (1 , valign = "top" , latex_hline = 'none' ) %>%
group_rows ("Primary outcome" , 1 , nrow (comp_table)) %>%
group_rows ("Mortality" , nrow (comp_table) + 1 , nrow (comp_table) + nrow (mort_table)) %>%
group_rows ("Vasopressor/inotropes" , 11 , 15 ) %>%
group_rows ("Ventilation" , 16 , 20 )
save_tex_table (out, fn = "outcomes/primary/composite-breakdown" )
```
```{r}
#| label: tbl-po-composite-combinations
#| code-summary: Summarise the component combinations
#| tbl-cap: |
#| Observed outcome combinations.
levs <- c ("None" , "Death + Vas./Ino. + Ventilation" , "Death + Vas./Ino." , "Death + Ventilation" , "Vas./Ino. + Ventilation" ,
"Death" , "Vas./Ino." , "Ventilation" , "Unknown" )
tab <- all_dat %>%
filter_acs_itt () %>%
group_by (CAssignment) %>%
summarise (
` Composite outcomes ` = case_when (
PO == 0 ~ 0 ,
is.na (PO) ~ 99 ,
D28_death == 1 & ANY_vasop == 1 & ANY_vent == 1 ~ 1 ,
D28_death == 1 & ANY_vasop == 1 & ANY_vent == 0 ~ 2 ,
D28_death == 1 & ANY_vasop == 0 & ANY_vent == 1 ~ 3 ,
D28_death == 0 & ANY_vasop == 1 & ANY_vent == 1 ~ 4 ,
D28_death == 1 ~ 5 ,
ANY_vasop == 1 ~ 6 ,
ANY_vent == 1 ~ 7 ,
)) %>%
dplyr:: count (CAssignment, ` Composite outcomes ` ) %>%
mutate (` Composite outcomes ` = factor (` Composite outcomes ` , levels = c (0 : 7 , 99 ), labels = levs)) %>%
spread (CAssignment, n, fill = 0 ) %>%
mutate (Overall = C1 + C2 + C3 + C4) %>%
mutate (across (C1: Overall, ~ sprintf ("%i (%.0f)" , .x, 100 * .x / sum (.x))))
colnames (tab)[- 1 ] <- tab_head
colnames (tab)[1 ] <- "Composite outcomes, n ( \\ %)"
kable (
tab ,
align = "lrrrrr" ,
escape = F) %>%
kable_styling (font_size = 11 )
```
```{r}
colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )
out <- kable (
tab,
format = "latex" ,
align = "lrrrrr" ,
booktabs = TRUE ,
linesep = "" ,
escape = F) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position" )
save_tex_table (out, fn = "outcomes/primary/composite-combinations" )
```
# Descriptive Analyses
## By Interim Analysis
```{r}
intdat <- all_dat %>%
filter_acs_itt () %>%
mutate (
interim = as.character (case_when (
RandDate < get_interim_dates ()$ meet_date[1 ] ~ 1 ,
RandDate < get_interim_dates ()$ meet_date[2 ] ~ 2 ,
RandDate < get_interim_dates ()$ meet_date[3 ] ~ 3 ,
RandDate < get_interim_dates ()$ meet_date[4 ] ~ 4 ,
TRUE ~ 4
)),
CAssignment = factor (
CAssignment, labels = c ("Low" , "Intermediate" , "Low \n with aspirin" , "Therapuetic" ))
)
intcounts <- intdat %>%
dplyr:: count (interim) %>%
mutate (label = paste0 ("Interim " , interim, " (n = " , n, ")" ))
trtcounts <- intdat %>% dplyr:: count (CAssignment)
ovrtab <- intdat %>%
group_by (interim = "Overall" , CAssignment) %>%
summarise (
Randomised =
sprintf ("%i" , n ()),
Known =
sprintf ("%i (%.0f)" , sum (! is.na (PO)), 100 * mean (! is.na (PO))),
` Met primary outcome ` =
sprintf ("%i (%.0f)" , sum (PO, na.rm = T), 100 * mean (PO, na.rm = T)),
` Death ` =
sprintf ("%i (%.0f)" , sum (D28_death, na.rm = T), 100 * mean (D28_death, na.rm = T)),
` Vasopressor/inotropic support ` =
sprintf ("%i (%.0f)" , sum (ANY_vasop, na.rm = T), 100 * mean (ANY_vasop, na.rm = T)),
` Intensive respiratory support ` =
sprintf ("%i (%.0f)" , sum (ANY_vent, na.rm = T), 100 * mean (ANY_vent, na.rm = T))
) %>%
pivot_longer (` Randomised ` : ` Intensive respiratory support ` , names_to = "Outcome" ) %>%
pivot_wider (names_from = CAssignment, values_from = value, values_fill = "-" )
inttab <- intdat %>%
group_by (interim, CAssignment) %>%
summarise (
Randomised =
sprintf ("%i" , n ()),
Known =
sprintf ("%i (%.0f)" , sum (! is.na (PO)), 100 * mean (! is.na (PO))),
` Met primary outcome ` =
sprintf ("%i (%.0f)" , sum (PO, na.rm = T), 100 * mean (PO, na.rm = T)),
` Death ` =
sprintf ("%i (%.0f)" , sum (D28_death, na.rm = T), 100 * mean (D28_death, na.rm = T)),
` Vasopressor/inotropic support ` =
sprintf ("%i (%.0f)" , sum (ANY_vasop, na.rm = T), 100 * mean (ANY_vasop, na.rm = T)),
` Intensive respiratory support ` =
sprintf ("%i (%.0f)" , sum (ANY_vent, na.rm = T), 100 * mean (ANY_vent, na.rm = T))
) %>%
pivot_longer (` Randomised ` : ` Intensive respiratory support ` , names_to = "Outcome" ) %>%
pivot_wider (names_from = CAssignment, values_from = value, values_fill = "-" )
tab <- bind_rows (ovrtab, inttab)
savetab <- kable (
tab[, - 1 ],
format = "latex" ,
align = "lrrrr" ,
booktabs = TRUE
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
group_rows (paste0 ("Overall (n = " , sum (intcounts$ n), ")" ), 1 , 6 ) %>%
group_rows (intcounts$ label[1 ], 7 , 12 ) %>%
group_rows (intcounts$ label[2 ], 13 , 18 ) %>%
group_rows (intcounts$ label[3 ], 19 , 24 ) %>%
group_rows (intcounts$ label[4 ], 25 , 30 )
save_tex_table (savetab, fn = "outcomes/primary/primary-composite-interims" )
kable (
tab[, - 1 ],
format = "html" ,
align = "lrrrr" ,
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
group_rows (paste0 ("Overall (n = " , sum (intcounts$ n), ")" ), 1 , 6 ) %>%
group_rows (intcounts$ label[1 ], 7 , 12 ) %>%
group_rows (intcounts$ label[2 ], 13 , 18 ) %>%
group_rows (intcounts$ label[3 ], 19 , 24 ) %>%
group_rows (intcounts$ label[4 ], 25 , 30 )
```
## Age
```{r}
#| label: fig-age-po
#| code-summary: Relationship age to outcome
#| fig-cap: |
#| Relationship (logistic regression linear in age)
#| between age at entry and the primary outcome.
agedat <- acs_itt_dat %>%
dplyr:: count (PO, AgeAtEntry) %>%
spread (PO, n, fill = 0 ) %>%
mutate (
n = ` 0 ` + ` 1 ` + ` <NA> ` ,
p = ` 1 ` / (` 1 ` + ` 0 ` ))
agemod <- glm (
cbind (` 1 ` , ` 0 ` ) ~ AgeAtEntry,
data = agedat,
family = binomial ())
agedat <- agedat %>%
mutate (
ypred = predict (agemod, newdata = agedat, type = "response" )
)
p1 <- ggplot (agedat, aes (AgeAtEntry, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" )
p2 <- ggplot (agedat, aes (AgeAtEntry, p)) +
geom_point () +
geom_vline (xintercept = 60 , linetype = 2 ) +
geom_line (aes (y = ypred)) +
labs (y = "Proportion \n with outcome" , x = "Age at entry" )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "primary" )
ggsave (file.path (path, "outcome-age.pdf" ), p, height = 2 , width = 6 )
p
```
```{r}
#| label: tbl-age
#| code-summary: Relationship age \u2265 60 to outcome
#| tbl-cap: |
#| Relationship between age >= 60 and the primary outcome.
acs_itt_dat %>%
dplyr:: count (` Age 60+ ` = agegte60, PO) %>%
group_by (` Age 60+ ` ) %>%
spread (PO, n, fill = 0 ) %>%
mutate (
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
) %>%
kable (digits = 2 ) %>%
kable_styling (bootstrap_options = "striped" )
```
## Country
```{r}
#| label: fig-country-po
#| code-summary: Relationship country to outcome
#| fig-cap: Primary outcome by country.
tdat <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" )), PO) %>%
group_by (Country) %>%
spread (PO, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
)
p1 <- ggplot (tdat, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (tdat, aes (Country, p_1)) +
geom_point () +
labs (
y = "Proportion meeting \n primary outcome" ,
x = "Country of enrolment" ) +
ylim (0 , NA )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "primary" )
ggsave (file.path (path, "outcome-country.pdf" ), p, height = 2 , width = 6 )
p
```
## Site
```{r}
#| label: fig-site-po
#| code-summary: Relationship site to outcome
#| fig-cap: Primary outcome by site within country.
tdat <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (
Country_lab = Country,
Site_lab = fct_infreq (Location),
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" )),
Site = PT_LocationName,
PO) %>%
group_by (Country, Site) %>%
spread (PO, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
) %>%
ungroup ()
p1 <- ggplot (tdat, aes (Site_lab, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (tdat, aes (Site_lab, p_1)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_point () +
labs (
y = "Proportion meeting \n primary outcome" ,
x = "Site of enrolment" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p <- p1 / p2
path <- file.path ("outputs" , "figures" , "outcomes" , "primary" )
ggsave (file.path (path, "outcome-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-cal-po
#| code-summary: Relationship calendar date to outcome
#| fig-cap: |
#| Relationship between calendar date and the primary outcome.
caldat <- acs_itt_dat %>%
dplyr:: count (PO, yr = year (RandDate), mth = month (RandDate)) %>%
spread (PO, n, fill = 0 ) %>%
mutate (p = ` 1 ` / (` 1 ` + ` 0 ` ),
n = ` 1 ` + ` 0 ` + ` <NA> ` )
p1 <- ggplot (caldat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_point () +
labs (
y = "Proportion meeting \n primary outcome" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 ) +
ylim (0 , 0.25 )
p2 <- ggplot (caldat, aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p <- p2 | p1
path <- file.path ("outputs" , "figures" , "outcomes" , "primary" )
ggsave (file.path (path, "outcome-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Days since symptom onset
```{r}
#| label: fig-dsfs-po
#| code-summary: Relationship days since sympton onset to outcome
#| fig-cap: |
#| Relationship (logistic regression linear)
#| between days since symptom onset at entry and the primary outcome.
dsfsdat <- acs_itt_dat %>%
dplyr:: count (PO, dsfs) %>%
spread (PO, n, fill = 0 ) %>%
mutate (n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p = ` 1 ` / (` 1 ` + ` 0 ` ))
dsfsmod <- glm (
cbind (` 1 ` , ` 0 ` ) ~ dsfs,
data = dsfsdat,
family = binomial ())
dsfsdat <- dsfsdat %>%
mutate (
ypred = predict (dsfsmod, newdata = dsfsdat, type = "response" )
)
p1 <- ggplot (dsfsdat, aes (dsfs, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 7.5 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Days since first symptoms" )
p2 <- ggplot (dsfsdat, aes (dsfs, p)) +
geom_point () +
geom_vline (xintercept = 7.5 , linetype = 2 ) +
geom_line (aes (y = ypred)) +
labs (y = "Proportion \n with outcome" ,
x = "Days since first symptoms" )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "primary" )
ggsave (file.path (path, "outcome-dsfs.pdf" ), p, height = 2 , width = 6 )
p
```
## Anti-coagulation Domain
```{r}
#| label: tbl-domainC-po
#| code-summary: Relationship anti-coagulation to outcome
#| tbl-cap: Primary outcome by anti-coagulation intervention.
make_po_table (
all_dat %>%
filter_acs_itt ()
)
save_tex_table (make_po_table (
all_dat %>%
filter_acs_itt (),
"latex" ),
"outcomes/primary/anticoagulation-summary" )
```
# Analyses
## ACS-ITT
In these analyses, the population is restricted to participants who were randomised to the anti-coagulation domain.
### Treatment only
```{r}
#| label: sample-trtonly-acs-itt-model
#| code-summary: Data and sampling treatment only model
res <- fit_primary_model (acs_itt_nona_dat, logistic_mod, seed = 17250 )
names (res$ drws$ OR) <- intervention_labels2 ()$ CAssignment[- (1 : 2 )]
```
```{r}
#| label: fig-or-densities-trtonly-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
plot_or_densities (res$ drws$ OR)
```
```{r}
#| label: tbl-or-densities-trtonly-acs-itt-model
#| tbl-cap: Posterior summaries for odds ratio contrasts.
odds_ratio_summary_table (res$ drws$ OR)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary ("beta" )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
### Treatment plus covariates (no site or epoch)
Under this model a component is added for country (region), $\mathbb{I}[ \text{age}\geq 60 ] $, and intervention specific ineligibility (C3 only)
::: column-margin
New terms:
$$
\begin{aligned}
\rho_1 &= 0 \\
\rho_r &\overset{\text{iid}}{\sim} \text{Normal}(0, 1), r=2,...,R,\\
\beta_{\text{age}\geq 60} &\sim \text{Normal}(0, 2.5^2) \\
\beta_{\text{inelgC3}} &\sim \text{Normal}(0, 10)
\end{aligned}
$$
:::
```{r}
#| label: sample-no-site-epoch-acs-itt-model
res <- fit_primary_model (
acs_itt_nona_dat,
logistic_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 71236 )
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: fig-or-densities-noepochsite-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
plot_or_densities (res$ drws$ OR)
```
```{r}
#| label: tbl-or-densities-noepochsite-acs-itt-model
#| tbl-cap: Posterior summaries for odds ratio contrasts.
odds_ratio_summary_table (res$ drws$ OR)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary ("beta" )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
### Treatment plus covariates plus site (no epoch)
Under this model a component is added for site nested within country.
::: column-margin
New terms
$$
\begin{aligned}
\xi_{rs} &\overset{\text{iid}}{\sim}\text{Normal}(0,\sigma^2_{r}), s=1,...,S_r, r=1,...,R \\
\sigma_r &\overset{\text{iid}}{\sim} \text{Half-}t(3, 1).
\end{aligned}
$$
:::
```{r}
#| label: model-noepoch-acs-itt-fit
#| code-summary: Data and sampling no site or epoch model
res <- fit_primary_model (
acs_itt_nona_dat,
logistic_site_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 52765 ,
adapt_delta = 0.95 )
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: fig-or-densities-noepoch-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
plot_or_densities (res$ drws$ OR)
```
```{r}
#| label: tbl-or-densities-noepoch-acs-itt-model
#| tbl-cap: Posterior summaries for odds ratio contrasts.
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
#| label: fig-site-effects-noepoch-acs-itt-model
#| fig-cap: Site effect odds ratios.
plot_site_terms (
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" , "gamma_site" , "tau_site" ))
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["tau_site" ])
```
:::
### Treatment plus covariates plus epoch (no site)
::: column-margin
Replace site terms with epoch terms
$$
\begin{aligned}
\tau_1 &= 0 \\
\tau_t &= \tau_{t-1} + \sigma_\tau\epsilon_t \\
\epsilon_t &\overset{\text{iid}}{\sim} \text{Normal}(0,1) \\
\sigma_\tau &\sim \text{Half-}t(3, 1)
\end{aligned}
$$
:::
```{r}
#| label: model-nosite-acs-itt-fit
#| code-summary: Data and sampling no site model
res <- fit_primary_model (
acs_itt_nona_dat,
logistic_epoch_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 76405 ,
adapt_delta = 0.98 )
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: fig-or-densities-nosite-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
plot_or_densities (res$ drws$ OR)
```
```{r}
#| label: tbl-or-densities-nosite-acs-itt-model
#| tbl-cap: Posterior summaries for odds ratio contrasts.
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
#| label: fig-epoch-effects-nosite-acs-itt-model
#| fig-cap: Epoch effect odds ratios.
plot_epoch_terms (
exp (res$ drws$ gamma_epoch)
)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" , "gamma_epoch" , "tau_epoch" ))
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_epoch" ])
mcmc_trace (res$ drws["tau_epoch" ])
```
:::
### Primary Model
```{r}
#| label: model-primary-sampling
#| code-summary: Data and sampling primary model
res <- fit_primary_model (
acs_itt_nona_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 32915 ,
adapt_delta = 0.99
)
res$ fit$ save_object (file.path ("outputs" , "models" , "primary" , "acs-itt-primary.rds" ))
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: odds-ratio-summary-table-primary-model
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table (res$ drws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-summary-table" )
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
#| label: tbl-decision-quantities
#| tbl-cap: Decision quantity summaries
save_tex_table (
decision_quantity_summary_table (c ("Low-dose" = rvar (1 ), res$ drws$ OR[1 : 3 ]), "latex" ),
"outcomes/primary/primary-model-acs-itt-decision-table" )
decision_quantity_summary_table (
c ("Low-dose" = rvar (1 ), res$ drws$ OR[1 : 3 ])
)
```
```{r}
#| label: fig-or-densities-acs-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
p <- plot_or_densities (res$ drws$ OR[1 : 3 ]) +
labs (x = "Odds ratio (log scale)" , y = "Comparison" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "primary-model-acs-itt-odds-ratio-densities.pdf" )
ggsave (pth, p, width = 6 , height = 2.5 )
p
```
```{r}
#| label: odds-ratio-summary-primary-model-acs-itt-epoch-site-terms
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "primary-model-acs-itt-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" , "gamma_epoch" , "tau_epoch" , "gamma_site" , "tau_site" ))
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["gamma_epoch" ])
mcmc_trace (res$ drws["tau_epoch" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["tau_site" ])
```
:::
#### Posterior Predictive Checks
```{r}
#| label: primary-acs-itt-ppc
y_ppc <- res$ drws$ y_ppc
ppc_dat <- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc))
grp_ppc <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
mean_y = mean (PO),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
}
plot_grp_ppc <- function (dat, lab = "" ) {
dat %>%
ggplot (aes (y = grp, xdist = rvar_mean_y_ppc)) +
stat_interval (size = 2 ) +
geom_point (aes (x = mean_y, y = 1 : nrow (dat)), colour = "red" , shape = 23 ) +
labs (
x = "Posterior predictive \n proportion" ,
y = lab)
}
ppc_C <- grp_ppc (sym ("CAssignment" ))
ppc_ctry <- grp_ppc (sym ("country" ))
ppc_epoch <- grp_ppc (sym ("epoch" ))
ppc_site <- ppc_dat %>%
group_by (Country = country, grp = site_raw) %>%
summarise (
mean_y = mean (PO),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
p1 <- plot_grp_ppc (ppc_C, "Anticoagulation \n intervention" ) + labs (x = "" )
p2 <- plot_grp_ppc (ppc_ctry, "Country" ) + labs (x = "" )
p3 <- plot_grp_ppc (ppc_epoch, "Epoch" ) + labs (x = "" )
p4 <- plot_grp_ppc (ppc_site %>% filter (Country == "IN" ), "Sites \n India" ) + labs (x = "" )
p5 <- plot_grp_ppc (ppc_site %>% filter (Country == "AU" ), "Sites \n Australia" ) + labs (x = "" )
p6 <- plot_grp_ppc (ppc_site %>% filter (Country == "NP" ), "Sites \n Nepal" )
p7 <- plot_grp_ppc (ppc_site %>% filter (Country == "NZ" ), "Sites \n NZ" )
p <- ((p3 | p1 / p2) /
( (p4 | p5) / (p6 | p7) +
plot_layout (heights = c (3 , 1 )) ) ) +
plot_layout (heights = c (1 , 1.5 ), guides = "collect" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" ,
"primary-model-acs-itt-ppc.pdf" )
ggsave (pth, p, width = 6 , height = 5.5 )
p
```
## Sensitivity Analyses
### Sens: Treatment Contrasts
```{r}
res <- fit_primary_model (
acs_itt_nona_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
ctr = contr.treatment,
seed = 32915 ,
adapt_delta = 0.99
)
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| tbl-cap: Posterior summaries for model parameters under treatment contrast priors.
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
save_tex_table (
odds_ratio_summary_table (res$ drws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-trt-ctr-summary-table" )
```
```{r}
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" ,
"primary-model-acs-itt-trt-ctr-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Sens: Weakly Informative Priors
As a prior sensitivity assessment, the primary model is re-fit under weakly informative priors (Normal(0, 10)) for all fixed-effect parameters.
```{r}
res <- fit_primary_model (
acs_itt_nona_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 10 , 10 , 10 ),
seed = 32915 ,
adapt_delta = 0.99
)
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| tbl-cap: Posterior summaries for model parameters under treatment contrast priors.
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
p
```
### Sens: FAS-ITT
In these analyses, the population is includes all participants who were randomised to the platform, even if not randomised into the anticoagulation domain.
```{r}
#| label: sample-primary-fas-itt-model
res <- fit_primary_model (
fas_itt_nona_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 32915 ,
adapt_delta = 0.99
)
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| label: fig-or-densities-primary-fas-itt-model
#| fig-cap: Posterior densities for odds ratio contrasts.
#| echo: false
plot_or_densities (res$ drws$ OR)
```
```{r}
#| label: tbl-or-densities-primary-fas-itt-model
#| tbl-cap: Posterior summaries for odds ratio contrasts.
#| echo: false
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
#| label: odds-ratio-summary-table-primary-model-fas-itt-save
save_tex_table (
odds_ratio_summary_table (res$ drws$ OR, "latex" ),
"outcomes/primary/primary-model-fas-itt-summary-table" )
```
```{r}
#| label: odds-ratio-summary-primary-model-fas-itt-epoch-site-terms
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" , "primary-model-fas-itt-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Sens: Worst-Case
Under this scenario, all participants with unknown primary outcome status were were to assumed to have met the primary outcome.
```{r}
make_po_table (acs_itt_wc_dat)
save_tex_table (
make_po_table (acs_itt_wc_dat, "latex" ),
"outcomes/primary/acs-itt-worst-case-anticoagulation-summary" )
```
```{r}
#| label: sample-primary-acs-itt-wc-model
res <- fit_primary_model (
acs_itt_wc_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 32915 ,
adapt_delta = 0.99
)
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| fig-cap: Posterior densities for odds ratio contrasts.
#| echo: false
plot_or_densities (res$ drws$ OR)
```
```{r}
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table (res$ drws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-worst-case-summary-table" )
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" )))
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" ,
"primary-model-acs-itt-worst-case-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Sens: Best-Case
Under this scenario, all participants with unknown primary outcome status were were to assumed to have **not** met the primary outcome.
```{r}
save_tex_table (
make_po_table (acs_itt_bc_dat, "latex" ),
"outcomes/primary/acs-itt-best-case-anticoagulation-summary" )
make_po_table (acs_itt_bc_dat)
```
```{r}
#| label: sample-primary-acs-itt-bc-model
res <- fit_primary_model (
acs_itt_bc_dat,
primary_mod,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
beta_sd_var = c (10 , 2.5 , 1 , 1 ),
seed = 32915 ,
adapt_delta = 0.99
)
names (res$ drws$ OR) <- c (
intervention_labels2 ()$ CAssignment[- (1 : 2 )],
"Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" )
```
```{r}
#| fig-cap: Posterior densities for odds ratio contrasts.
#| echo: false
plot_or_densities (res$ drws$ OR)
```
```{r}
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table (res$ drws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-best-case-summary-table" )
odds_ratio_summary_table (res$ drws$ OR)
```
```{r}
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (
res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" )))
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" ,
"primary-model-acs-itt-best-case-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
### Sens: Concurrent Randomisation
#### Intermediate-dose
```{r}
#| label: tbl-acs-itt-c2-summary-table
#| tbl-cap: Primary outcome summary for participants randomised to either
#| low-dose or intermediate dose arms (ACS-ITT-intermediate).
save_tex_table (
make_po_table (
acs_itt_concurc2_dat,
"latex" ),
"outcomes/primary/acs-itt-c2-summary" )
make_po_table (acs_itt_concurc2_dat)
```
```{r}
#| label: fit-model-acs-itt-c2
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc2_nona_dat,
contrasts.arg = list (CAssignment = ctr))
y <- acs_itt_concurc2_nona_dat$ PO
sdat <- list (
N = nrow (X),
K = ncol (X),
X = X,
y = y,
beta_sd = c (2.5 , 1 , 2.5 , 1 , 1 )
)
snk <- capture.output (
sfit <- logistic_mod$ sample (
data = sdat,
seed = 3274 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500
)
)
sdrws <- as_draws_rvars (sfit$ draws ("beta" ))
names (sdrws$ beta) <- colnames (X)
sdrws$ beta_trt <- as.vector (
transform_coding (cbind (1 , ctr (2 )), cbind (1 , contr.treatment (2 ))) %**% sdrws$ beta[1 : 2 ]
)
sdrws$ OR <- c (
"Intermediate-dose" = exp (sdrws$ beta_trt[2 ]),
"Age \u2265 60" = exp (sdrws$ beta[3 ]),
"Australia/New Zealand" = exp (sdrws$ beta[4 ]),
"Nepal" = exp (sdrws$ beta[5 ])
)
```
```{r}
#| label: tbl-acs-itt-c2-primary-model-summary
#| tbl-cap: Parameter posterior summaries.
save_tex_table (
odds_ratio_summary_table (sdrws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-c2-summary-table" )
odds_ratio_summary_table (sdrws$ OR)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
res$ fit$ summary (c ("beta" ))
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
```
:::
#### Low-dose with aspirin
```{r}
#| label: tbl-acs-itt-c3-summary-table
#| tbl-cap: Primary outcome summary for participants randomised whilst the
#| low-dose plus aspirin arm was available (ACS-ITT-aspirin).
save_tex_table (
make_po_table (
acs_itt_concurc3_dat,
"latex" ),
"outcomes/primary/acs-itt-c3-summary" )
make_po_table (acs_itt_concurc3_dat)
```
```{r}
#| label: fit-model-acs-itt-c3
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc3_nona_dat,
contrasts.arg = list (CAssignment = contr.orthonorm))
y <- acs_itt_concurc3_nona_dat$ PO
sdat <- list (
N = nrow (X),
K = ncol (X),
X = X,
y = y,
beta_sd = c (2.5 , 1 , 1 , 2.5 , 1 )
)
snk <- capture.output (
sfit <- logistic_mod$ sample (
data = sdat,
seed = 3274 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500
)
)
sdrws <- as_draws_rvars (sfit$ draws ("beta" ))
names (sdrws$ beta) <- colnames (X)
sdrws$ beta_trt <- as.vector (
transform_coding (cbind (1 , contr.orthonorm (3 )), cbind (1 , contr.treatment (3 ))) %**% sdrws$ beta[1 : 3 ]
)
sdrws$ OR <- c (
"Intermediate-dose" = exp (sdrws$ beta_trt[2 ]),
"Low-dose with aspirin" = exp (sdrws$ beta_trt[3 ]),
"Age \u2265 60" = exp (sdrws$ beta[4 ]),
"Australia/New Zealand" = exp (sdrws$ beta[5 ])
)
```
```{r}
#| label: tbl-acs-itt-c3-primary-model-summary
#| tbl-cap: Parameter posterior summaries.
save_tex_table (
odds_ratio_summary_table (sdrws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-c3-summary-table" )
odds_ratio_summary_table (sdrws$ OR)
```
#### Therapeutic-dose
```{r}
#| label: tbl-acs-itt-c4-summary-table
#| tbl-cap: Primary outcome summary for participants randomised under protocol 5.0 (ACS-ITT-therapeutic).
save_tex_table (
make_po_table (
acs_itt_concurc4_dat,
"latex" ),
"outcomes/primary/acs-itt-c4-summary" )
make_po_table (acs_itt_concurc4_dat)
```
```{r}
#| label: fit-model-acs-itt-c4
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc4_nona_dat,
contrasts.arg = list (CAssignment = contr.orthonorm))
y <- acs_itt_concurc4_nona_dat$ PO
sdat <- list (
N = nrow (X),
K = ncol (X),
X = X,
y = y,
beta_sd = c (2.5 , 1 , 1 , 2.5 , 1 , 1 )
)
snk <- capture.output (
sfit <- logistic_mod$ sample (
data = sdat,
seed = 3274 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500
)
)
sdrws <- as_draws_rvars (sfit$ draws ("beta" ))
names (sdrws$ beta) <- colnames (X)
sdrws$ beta_trt <- as.vector (
transform_coding (cbind (1 , contr.orthonorm (3 )), cbind (1 , contr.treatment (3 ))) %**% sdrws$ beta[1 : 3 ]
)
sdrws$ OR <- c (
"Intermediate dose" = exp (sdrws$ beta_trt[2 ]),
"Therapeutic dose" = exp (sdrws$ beta_trt[3 ]),
"Age \u2265 60" = exp (sdrws$ beta[4 ]),
"India" = exp (sdrws$ beta[5 ]),
"Australia/New Zealand" = exp (sdrws$ beta[6 ])
)
```
```{r}
#| label: tbl-acs-itt-c4-primary-model-summary
#| tbl-cap: Parameter posterior summaries.
save_tex_table (
odds_ratio_summary_table (sdrws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-itt-c4-summary-table" )
odds_ratio_summary_table (sdrws$ OR)
```
## Per-Protocol Analysis
In the per-protocol analysis, participants who did not satisfy the intervention protocol (as determined by blinded review) are excluded from the analysis set.
```{r}
all_pp %>%
filter (ENR_rec == 1 ) %>%
dplyr:: count (Reason)
```
```{r}
acs_itt_dat %>%
left_join (pp, by = "StudyPatientID" ) %>%
dplyr:: count (CAssignment, Reason) %>%
spread (CAssignment, n, fill = 0 )
acs_itt_dat %>%
left_join (pp, by = "StudyPatientID" ) %>%
dplyr:: count (country, Reason) %>%
spread (country, n, fill = 0 )
```
```{r}
#| label: tbl-acs-pp-summary-table
#| tbl-cap: Primary outcome summary for ACS-PP set.
save_tex_table (
make_po_table (
acs_pp_dat,
"latex" ),
"outcomes/primary/acs-pp-outcome-summary" )
make_po_table (acs_pp_dat)
```
```{r}
#| label: model-primary-sampling-acs-pp
#| code-summary: Data and sampling primary model
mdat <- make_primary_model_data (
acs_pp_nona_dat,
c ("inelgc3" , "agegte60" , "ctry" ),
c (10 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- primary_mod$ sample (
data = mdat,
seed = 81241 ,
refresh = 0 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
chains = 8 ,
parallel_chains = 8 ,
adapt_delta = 0.99
))
mdrws <- as_draws_rvars (mfit$ draws (c (
"beta" ,
"gamma_site" , "tau_site" ,
"gamma_epoch" , "tau_epoch" , "y_ppc" )))
names (mdrws$ beta) <- colnames (mdat$ X)
# Transformed samples
epoch_map <- acs_pp_nona_dat %>% dplyr:: count (epoch, epoch_lab)
site_map <- acs_pp_nona_dat %>% dplyr:: count (site_num, site)
site_facet <- factor (mdat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
Ca <- mdat$ ctrA
Cc <- mdat$ ctrC
mdrws$ Acon <- Ca %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- Cc %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]
mdrws$ Atrt <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ OR <- setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" ))
mdrws$ OR <- c (
mdrws$ OR,
"Ineligible aspirin" = exp (mdrws$ beta[7 ]),
"Age \u2265 60" = exp (mdrws$ beta[8 ]),
setNames (exp (mdrws$ beta[9 : 10 ]), c ("Australia/New Zealand" , "Nepal" )))
mdrws$ ORepoch <- setNames (exp (mdrws$ gamma_epoch), epoch_map$ epoch_lab)
mdrws$ ORsites <- setNames (exp (mdrws$ gamma_site), site_map$ site)
```
```{r}
#| label: tbl-odds-ratio-summary-table-primary-model-acs-pp
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table (mdrws$ OR, "latex" ),
"outcomes/primary/primary-model-acs-pp-summary-table" )
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
p <- plot_epoch_site_terms (mdrws$ ORepoch, mdrws$ ORsites, site_facet)
pth <- file.path ("outputs" , "figures" , "outcomes" , "primary" ,
"primary-model-acs-pp-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```